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We consider the coaction of two distinct noise sources on the activation process of a single and two interacting 
excitable units, which are mathematically described by the Fitzhugh-Nagumo equations. We determine the most 
probable activation paths around which the corresponding stochastic trajectories are clustered. The key point 
lies in introducing appropriate boundary conditions that are relevant for a class II excitable unit, which can be 
immediately generalized also to scenarios involving two coupled units. We analyze the effects of the two noise 
sources on the statistical features of the activation process, in particular demonstrating how these are modified 
due to the linear/nonlinear form of interactions. Universal properties of the activation process are qualitatively 
discussed in the light of a stochastic bifurcation that underlies the transition from a stochastically stable fixed 
point to continuous oscillations. 


I. INTRODUCTION 

Excitability is a dynamical feature shared by nonlinear sys¬ 
tems that lie in vicinity of bifurcation underlying transition 
from stationary state toward the sustained periodic activity (U] . 
Study of models where stochastic excitable dynamics is cru¬ 
cial for shaping local or collective behavior is a rapidly devel¬ 
oping field, gaining relevance in terms of theory as well as for 
offering qualitative insights into a variety of biological M 
and inorganic systems |0-lal • 

Excitability is manifested in the existence of a narrow range 
of stimuli magnitudes where marginally different perturba¬ 
tions may cause the system to generate two qualitatively dis¬ 
tinct types of responses. While smaller perturbations give 
rise to small-amplitude (linear) responses, the slightly larger 
perturbations may elicit pulse-like, large-amplitude excitation 
loops. The latter are comprised of activation and relaxation 
stage, whereby the spike profile remains independent on the 
form of applied perturbation. If one interprets perturbation in 
terms of setting the particular initial conditions, then the ex¬ 
citability feature is reflected in the point that the system shows 
strong sensitivity to initial conditions within a small domain 
of relevant values. Consequently, behavior of excitable sys¬ 
tems is heavily susceptible to noise l|3l, whose influence may 
at least in part be understood as excitability amplification. 

Given the strong likelihood of encountering systems which 
are influenced by combined action of variations in the envi¬ 
ronment and the fluctuations of the internal parameters, it is 
justified to analyze models where different sources of noise 
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affect the dynamics of multiple variables. This point, together 
with the fact that the existence of pulse-like excitations re¬ 
quires a flow with dimension no less than two ([^, suggests 
that the system with two variables, each subjected to stochas¬ 
tic perturbation, may be considered a paradigmatic setup. Eor 
example, two sharply separated characteristic time scales and 
two sources of noise acting independently on the fast and the 
slow variable are relevant ingredients for the description of 
a typical neuron HEl or laser dynamics i]. Eor neurons, 
stochastic term affecting the fast variable may account for 
synaptic noise due to random arrival of spikes from a large 
number of afferents, whereas the stochastic component in the 
slow variable dynamics may be seen as internal noise, asso¬ 
ciated to thermal fluctuations in the opening of the ion-gating 
channels. Consistent with this, we refer to stochastic term 
added to the dynamics of the fast (slow) variable as external 
(internal) noise. 

Erom the theoretical viewpoint, a comprehensive insight 
into the noise-driven activation process is fundamental for un¬ 
derstanding excitability. The problem we focus on concerns 
the activation process in a single or two coupled excitable 
Eitzhugh-Nagumo {FHN) elements driven by two indepen¬ 
dent sources of noise. The main reason for considering the 
FHN model lies in its generic character, viz. the fact that 
it is canonical for class II excitable systems, which involve an 
almost continuous transition between the small- and the large- 
amplitude excitations. Eor conceptual reasons and in view of 
potential applications, only the spiking responses are relevant 
and of interest to us. It follows that the formulation of activa¬ 
tion event has to be adapted to class II excitable systems, such 
that it satisfies two criteria; (i) it should warrant a clear-cut 
distinction of the spiking response from the small-amplitude 
excitation and (ii) it should allow an immediate generaliza¬ 
tion from the case of a single unit to that of two interacting 
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units. We stress that in order to meet both criteria, it is crucial 
to introduce an appropriate terminating boundary set for the 
activation problem, as demonstrated in the paper. 

Note that the term activation is inherited from, but it is not 
used in the same sense as in a typical escape problem 10 [ni- 
Q. To explain the differences, one first invokes a general 
remark that for excitable systems, two possible types of re¬ 
sponse derive from the proximity to bifurcation rather than 
the bistable dynamics, so that the threshold-like behavior is 
not associated to a genuine separatrix (saddle structure), but 
rests on the coaction of nonlinearity and the sharp separation 
between the system’s characteristic time scales. Thus, resolv¬ 
ing the origin and character of threshold-like behavior for ex¬ 
citable elements may be linked to an unstable fixed point or, 
as in case of a FHN unit, to a quasi-separatrix (’’ghost sep¬ 
aratrix”) ifT^ . This is apparently different from the typical 
escape problem scenario. 

Another important point is that we tie the activation ex¬ 
clusively to spiking response, so that the event where the 
phase point reaches the quasi-separatrix per se is insufficient 
to count as activation. In other words, crossing the quasi- 
separatrix does not present a discriminative condition between 
the small- and the large-amplitude excitations. This is why we 
introduce novel terminating boundary conditions, rather than 
just consider an extension of the escape problem to excitable 
systems s. One should also point out that our approach is 
conceptually different from the earlier work concerning noisy 
excitable neurons where the terminating boundary 

does not constitute a set, but rather a unique boundary point, 
introduced as an arbitrary threshold independent on the struc¬ 
ture of phase space. 

We analyze the activation problem from two angles, one 
by determining the most probable activation paths (MPAPs), 
and the other by examining the statistical features of the ac¬ 
tivation process. An important result consists in determining 
the MPAPs for a single and two coupled units under dif¬ 
ferent ratios of external vs internal noise Di/02- In case of 
two units, the topology of the obtained trajectories is shown 
to qualitatively depend on the form of coupling. The statis¬ 
tics of a single- and two-unit activation process is character¬ 
ized by examining how the time-to-first-pulse (T FP) t aver¬ 
aged over different stochastic realizations and the associated 
coefficient of variation R depend on noise. Both t{Di,D 2 ) 
and R{Di,D 2 ) are found to display universal behavior for all 
considered setups, whereby the transition between two of the 
observed r regimes is associated to the fact that a unit under¬ 
goes stochastic bifurcation induced by Di and D 2 . For two 
units, the form of coupling is shown to have a nontrivial effect 
on correlation of the individual mean TFPs, which we relate 
to synchronization properties of the time series for the given 
parameter set. 

The paper is organized as follows. In Sec. |II] we present 
the details of the model, focusing on the threshold-like be¬ 
havior of a single unit and the results of bifurcation analysis 
for the coupled units. Section Hill concerns the case of a sin¬ 
gle unit subjected to external and internal noise. Having laid 
out the details of the method used to determine the MPAPs, 
cf. subsection llll A1 we examine how the topological features 


of the MPAPs depend on the pertaining noise intensities. 
Apart from relating the properties of t{D 1 , 02 ) dependence 
to the onset of stochastic bifurcation from the stochastically 
stable fixed point to the stochastically stable limit cycle, we 
also introduce an approximation to explicitly demonstrate that 
Di and D 2 make substantially different impact on the mean 
TFPs. Section HVl contains the results for two units inter¬ 
acting via the linear or nonlinear couplings. It is analyzed 
how the different form of coupling affects the profile of the 
respective MPAPs, the stochastic bifurcation as well as the 
correlation of single unit mean TFPs. Section IVl provides a 
summary of the main results. 

II. BACKGROUND ON THE APPLIED MODEL 
A. Dynamics of a single excitable unit 

As a paradigm for excitable systems, we consider the 
FHN model, so that the dynamics of a single unit is given 
by 

dx = fx{x, y) = [x — x^/i — y]dt + y/2DidWi 

dy = fy{x,y) = e{x + b)dt+y^ 2 D^dW 2 . (1) 

The model is canonical for type II excitability class, meaning 
that the equilibrium lies in vicinity of the direct supercritical 
Hopf bifurcation. The latter is controlled by the excitability 
parameter b, whereby the critical value is |6| = 1. For |&| > 1, 
the system possesses a unique stable equilibrium (Xeq, yeq) = 
(—6, —b + b^/3), whereas for |6| < 1 the oscillatory state sets 
in. Given the symmetry of the system O, the analysis may be 
confined to case 6 > 0 without loss of generality. The unit in 
subcritical state displays excitable behavior if b is kept close 
to bifurcation threshold. In this paper, we fix 6 = 1.05. 

The other important ingredient of the model is the sharp 
separation between the characteristic time scales of the acti¬ 
vator and the recovery variable. The fast-slow dynamics is 
facilitated by setting e to a small value (e = 0.05). So far, 
FHN model has been applied in describing the ^namics 
of electrochemical reactions Q and cardiac cells 0, but is 
best known for its role in the field of neuroscience CHI]. In 
the latter context, the fast variable may be viewed as analo¬ 
gous to the neuron membrane potential, whereas the action 
of the slow variable may qualitatively be compared to that of 
K+ ion-gating channels yj]. Regarding the impact of ran¬ 
dom perturbations, we are interested into how the activation 
of units is shaped by two independent sources of noise. In ([1]), 
the stochastic effects are represented by the Wiener processes 
whose increments satisfy (dWi) = 0 and (dWidWj) = dtSij 
for i,j G 1,2. The dynamics of an excitable unit may be sum¬ 
marized as follows. In absence of perturbation, the selected 
parameter values are such that the system lies at equilibrium. 
Kicked by the perturbation, the unit may either display a small 
amplitude response, whereby the phase point rapidly decays 
back to equilibrium, or may exhibit a large excursion, settling 
to equilibrium only after the phase point has traversed the or¬ 
bit corresponding to a complete oscillation cycle. 
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FIG. 1. (Color online) Phase plane analysis for a FHN excitable 
unit. Equilibrium (EQ) lies at the intersection of the nullclines 
fx{x,y) = 0 and fy{x) = 0. The x-nullcline comprises three 
branches. In the singular limit e ^ 0, the spiking branch Ss and the 
refractory branch Sr are attractive, whereas the orbit separating the 
initial conditions that lead to small- or large-amplitude excitations, 
SAE and LAE respectively, contains the middle unstable branch. 
For finite e, the boundary between the two sets of initial conditions 
foliates into a thin layer of canard-like trajectories (CNRD) referred 
to as ’’ghost separatrix”. Illustrated trajectories are obtained by fixing 
the X initial condition to xo = —3, whereas yo is varied. Extreme 
sensitivity to initial conditions in vicinity of ghost separatrix is cor¬ 
roborated by the fact that a difference in yo of the order of 10” 
gives rise to a different form of excitation. The parameters of FHN 
model are fixed to e = 0.05, b = 1.05. 


For the proper statement of our problem, it is important 
to first consider the geometric interpretation of the unit’s dy¬ 
namics. Two particular issues are addressed, one related to 
the existence and the structure of the boundary between the 
initial conditions resulting in small- or large-amplitude re¬ 
sponses, whereas the other concerns the distinction between 
the roles of Di and D 2 in triggering the pulse emission. Re¬ 
garding the first point, we summarize the results on the de¬ 
terministic (noiseless) version of the system O obtained by 
the method of phase plane analysis, which rests on the sin¬ 
gular perturbation theory. In the limit e 0 , O turns into 
a one-dimensional system x = fx{x, y) under the constraint 
y — 0, meaning that y may be viewed as a fixed parameter. 
For small but finite e, x quickly relaxes to the value given by 
fx{x, y) = 0 , the condition outlining the curve referred to as 
the X (cubic) nullcline or the slow manifold, see Fig. [T] The x 
nullcline is comprised of three branches, whose stability fea¬ 
tures derive from the singular limit e —0. While the refrac¬ 
tory Sr and the spiking branch Ss may then be regarded as 
attractors, the middle branch is unstable and is a part of sep¬ 
aratrix between the attractive branches. The key point is that 
for small but finite e, the structure of the boundary and the 
related threshold behavior is mostly inherited from the singu¬ 
lar limit. The distinction is that finite e induces foliation of 
the boundary around the maximum of the x-nullcline. What 


has explicitly been shown by the so-called ’’blow-up method” 
ll^ . is that the boundary between the initial conditions lead¬ 
ing to Ss or Sr is not given by a single line, but rather by a 
thin layer made up of an infinite family of canard-like trajecto¬ 
ries of system ([T]l. The latter further implies that the boundary 
constitutes an invariant set. Boundary layer may still be per¬ 
ceived as a single line, a kind of ’’ghost separatrix’ ’ in, be¬ 
cause at distances d >> e from the fold point (1, 2/3), all the 
constituent trajectories become virtually indistinguishable. 

As for the qualitative interpretation of the effects of noise, 
it is evident that the noise term added to the slow-variable 
dynamics may shift the position of the y-nullcline. In other 
words, internal noise is capable of translating the fixed point 
from the stable refractory to the unstable (middle) branch of 
the X— nullcline, which temporarily switches the system from 
excitable to oscillatory state. While the impact of D 2 can be 
understood by geometric analysis, there is no analogous inter¬ 
pretation in case of Di. 

B. Dynamics of a couple of excitable units 

In case of a pair of coupled units, we consider two distinct 
setups, one where the units interact via the symmetrical lin¬ 
ear couplings, and the other involving couplings given by the 
nonlinear threshold-like function. 

The dynamics of a couple of FHN units interacting via 
linear couplings is given by 

dxi = [xi — a;f/3 — yi]dt + \/2DidWl + c[xi — Xj\dt 

dyi = €{xi + b)dt -F (2) 

where i,j G {l,2},i 7 ^ j specify the individual units. The 
random perturbations are such that the terms acting on differ¬ 
ent elements are supposed to be uncorrelated (dlF/dlT/) = 
0, fc,( G 1,2. Regarding the system parameters, note that 
the units are assumed to be identical (same h and e), while 
the symmetrical couplings are characterized by the coupling 
strength c. To see how the unit’s excitability feature is mod¬ 
ified in presence of interaction, we make a brief summary 
of the results of the bifurcation analysis carried out on the 
noise-free version of the system (O. The first remark is that 
the equilibrium is located at {xi,yi,X 2 ,y 2 ) = {—b,—b + 
6^/3, —b, —b + &^/3), whereby its stability is determined by 
the two pairs of complex conjugate characteristic exponents 
which satisfy ^ 1^2 = [1 — + 2c± ^(1 — b^ + 2c)^ — 4e]/2 

and = [1 — 6 ^ ± y/(l — 6 ^)^ — 4e]/2. Given that the sin¬ 
gle unit parameters are kept fixed, the local stability of equi¬ 
librium is changed under variation of the coupling strength. 
In particular, it may be demonstrated that the system under¬ 
goes a direct supercritical Hopf bifurcation at the critical value 
ch.l = (^^~l)/2. At this point, the stability of equilibrium is 
lost and the units are no longer in the excitable regime. How¬ 
ever, we stress that the complete picture on the system dy¬ 
namics cannot be gained from the local bifurcation analysis 
alone, because even before the excitability feature is lost (for 
c = cr^r), it is modified due to a global fold-cycle bifurca¬ 
tion which occurs at crc < cr.l- Thus, under increasing c 
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the system exhibits three types of characteristic behavior: (i) 
the ’’proper” excitable regime for c < cpc, {H) the regime of 
’’generalized excitability” for c G {cfCtCh,l) and (in) the 
oscillatory state at c > ch.l- The case (ii) involves coexis¬ 
tence between the fixed point and the large limit cycle, whose 
basins of attraction are separated by the stable manifold of the 
saddle cycle. In strict terms, such a scenario does not con¬ 
form to Izhikevich’s definition of excitability, because instead 
of relaxing to equilibrium after the large excursion, the sys¬ 
tem displays continuous oscillations. Nevertheless, some au¬ 
thors consider such a behavior as excitable or excitable-like, 
and one should also note that the oscillation orbit may still 
pass quite close to equilibrium, depending on the position, 
size and the manifold structure of the saddle cycle. Note that 
above ch.l, the large limit cycle created in the global event 
survives as the only attractor, because the incipient cycle born 
via Hopf bifurcation grows only until colliding with the pre¬ 
existing saddle-cycle, such that the two get annihilated in the 
inverse fold-cycle bifurcation. In other words, the properties 
of the limit cycle attractor both below and above ch,l are de¬ 
termined by the global bifurcation. 

Having summarized the results of the bifurcation analysis 
for the setup involving two units coupled via linear function, 
we turn to the scenario where the units subjected to external 
and internal noise interact in a nonlinear fashion. The dynam¬ 
ics of the units then obeys 

dxi = [xi — xf/3 — yi]dt + \J 2DidWl + caictan{xj + b)dt 
dyi = e(xi -f b)dt + ^/w^dW^, (3) 

The interaction terms have a form of a threshold-like func¬ 
tion, whose argument is defined relative to the correspond¬ 
ing unit’s equilibrium Xj — xj^eq = xj + b. This way, 
the impact of the state Xj is felt stronger if it lies further 
away from the equilibrium. The stability of equilibrium is 
determined by the four roots of the characteristic equation, 
which appear as two pairs of complex conjugates p-i 2 = 

[1 - b^ + c± ^{l-b^ + c)^-Ae]/2, /i 3.4 = [1 - b^ - 
c± y'jl — 6 ^ — c)2 — 4e]/2. Again, it may be shown that the 
system (O undergoes a direct supercritical Hopf bifurcation 
at ch.nl = 6 ^ — 1. Note that we confine the analysis to the 
case c > 0. Unlike the setup based on the linear coupling, 
here one does not encounter the global bifurcation controlled 
by c. In other words, the equilibrium is stable and the units 
are in excitable regime for c < ch,nl, whereas the system is 
in oscillatory state for c > ch,nl- 


III. ACTIVATION PROCESS IN AN EXCITABLE UNIT 
DRIVEN BY EXTERNAL AND INTERNAL NOISE 

In the following subsection we introduce the numerical 
method applied to determine the MPAPs. The method is 
illustrated in case of a single FHN unit, but it can readily 
carry over to the process of hrst pulse emission for two in¬ 
teracting units. In subsections IIIIBl and IIV Al the topolog¬ 
ical features of the pertaining trajectories will be analyzed 


in reference to stochastic bifurcation, viz. for noise intensi¬ 
ties substantially below, near or above the critical domain of 
(Di, D 2 ) values. While the main focus lies with the MPAPs 
explicitly obtained from stochastic simulations, the extended 
discussion will also concern the trajectories generated as so¬ 
lutions of the effective Hamiltonian equations under boundary 
conditions relevant for the process of first pulse emission. We 
stress that the use of such equations is distinct from the stan¬ 
dard Hamiltonian approach which yields optimal trajectories 
in a typical escape problem, or the generalization of an escape 
problem to an excitable FHN unit. In order to set up the 
discussion carried out in subsections IIII B I and IIV Al the fol¬ 
lowing subsection addresses the precise role of these effective 
equations and the ensuing trajectories. 


A. Method applied to determine the MPAPs 

The problem of obtaining the MPAPs for excitable units 
in general comprises two issues. One issue concerns how the 
terminating boundary conditions are specified, whereas the 
other relates to the particular details of the method. A com¬ 
mon ingredient in previous approaches to analysis of activa¬ 
tion process in excitable units has been to draw an analogy to 
motion of a particle in a one-dimensional potential perturbed 
by noise. Reduction to a one-dimensional problem naturally 
utilizes the decomposition of the system dynamics to fast and 
slow motions. Though such approximate methods are not in¬ 
tended to trace the unit’s most likely activation paths, the use 
of Fokker-Planck formalism still allows one to gain insight 
into certain statistical features of the activation process, in¬ 
cluding the mean activation time and its variance iiniiii. 
Nevertheless, an inherent drawback consists in the lack of 
ability to account for the simultaneous influence of perturba¬ 
tions added to both the slow and the fast component. In con¬ 
ceptual terms, the key problem lies in the fashion by which 
the escape from the stationary state is precisely defined. In 
particular, instead of associating the terminating boundary to 
the structure of the system’s phase space, the activation event 
is considered as a crossing of a predefined threshold, typically 
coinciding with the fold point at the minimum of the slow 
manifold. Compared to such methods, the approach we apply 
is preferred because ii) it introduces a unique definition of 
the activation path consistent with the structure of the phase 
space, (ii) one may consider the coaction of random pertur¬ 
bations added to both the fast and the slow subsystem, and 
(in), one is able to explicitly determine the MPAPs around 
which the different stochastic realizations are clustered. 

Before proceeding to the details of the numerical method 
implemented for calculation of the MPAPs, let us specify 
the boundary conditions relevant for the problem of first pulse 
emission in case of an excitable FHN unit. In particular, for 
the noise-driven system O, we consider the stochastic fluc¬ 
tuation paths in configuration space (x, y) that emanate from 
the deterministic fixed point {Xeq,yeq) = (- 6 ,-5 -f 6^/3) 
and terminate at the spiking branch of the cubic ((x)) null- 
cline. Given selection of the terminating boundary set derives 
from the fact that the spiking branch of the cubic nullcline 
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defines the spike profile for the deterministic limit cycle in 
the superthreshold regime b < 1, while the analogous point 
holds for the noise-induced oscillations in the excitable regime 
& > 1. In these terms, it has been verified that the profiles of 
spikes for an arbitrary combination of relevant {Di, D 2 ) val¬ 
ues can be approximated with sufficient accuracy by the orbit 
corresponding to the deterministic limit cycle characterizing 
the supercritical state. Note that the terminating boundary set 
we introduce is distinct from the one in ifisll . where a con¬ 
ceptually distinct problem has been considered. In particular, 
the aim in that paper has been to extend the standard escape 
problem to the case of an excitable FHN unit, which has 
been achieved by fixing the ’’ghost separatrix” as the termi¬ 
nating boundary set. The advantage of the boundary condi¬ 
tions we have adopted is that they can readily be generalized 
to the case of two interacting excitable units. In another pa¬ 
per, we shall further demonstrate that the relevance of a proper 
problem formulation which warrants that the small-amplitude 
response and the spiking response are clearly distinguished, 
as well as the related selection of the appropriate terminating 
boundary conditions, is even more pronounced when analyz¬ 
ing the noise-driven first pulse emission process for an assem¬ 
bly of excitable units. 

The details of the numerical method applied to determine 
the MPAPs are as follows. For the given {Di,D 2 ), we 
consider an ensemble of fluctuation paths {x{t),y{t)) which 
start from the deterministic fixed point {xeq, Ueq) at moment 
ti and satisfy the above stated terminating boundary condi¬ 
tions. The terminating time tf, as well as the associated coor¬ 
dinates {x{tf),y{tf)) are left unspecified. For the described 
ensemble, we study the statistics of the {x{t),y{t)) position 
of trajectories as a function of time ti <t <tf preceding the 
arrival to the terminating boundary set. In general, the idea is 
to sample the different stochastic realizations of activation tra¬ 
jectories in order to determine histograms of the path history 
reaching the terminating boundary set. Naturally, the recorded 
trajectories are characterized by the different tf times. Thus, 
the proper approach to characterize the statistics of the paths 
in confi gura tion space is to consider the prehistory probability 
density IlM , which concerns the distribution of paths ending 
at the specified boundary set. To obtain the former, one effec¬ 
tively sets the time when each stochastic realization terminates 
to f = 0, such that the behavior of the process during the ini¬ 
tiation of the pulse is observed by looking backward in time. 
This approach has been introduced in and has been ap¬ 
plied a number of times since ll^ . The prehistory probability 
distribution is defined as 

H{x,y,t)dxdy = Pr[x{t) G (x,x + dx),y(t) G {y,y + dy)\ 

— ^eq-;y{ii) — Veq]: 

ti <t < tf,x < Xb{tf),y < yb{tf). (4) 

The most probable path for the first pulse emission process 
is determined by collecting the points {xm{t), ym{t)) which 
correspond to the maximum of H{x,y,t) at any given mo¬ 
ment t. In other words, the MPAP for the given {Di, D 2 ) 
by definition coincides with the peak of as a function of 
time. In practice, the numerical method used to determine the 


prehistory probability density has involved dividing the (cc, y) 
phase space into a grid of 70 by 70 cells of length Ax = 0.048 
and width Ay = 0.09. Throughout the paper, the numerical 
integration of the appropriate set of stochastic equations is car¬ 
ried out by the Heun algorithm with the time step 5t = 0.002, 
whereas the averaging is performed over an ensemble of 5000 
different stochastic realizations of the first pulse emission pro¬ 
cess. 

In order to provide qualitative guidelines for an extended 
discussion, the MPAPs determined via the above described 
method for different setups and the different domains of 
noise intensities will be compared to trajectories obtained 
by integrating the set of effective Hamiltonian equations un¬ 
der boundary conditions relevant for the problem of first 
pulse emission. One recalls that in case of a typical escape 
problem, the Hamiltonian theory formulated in the extended 
variable-momentum state space may be used to explicitly ob¬ 
tain the optimal trajectories which coincide with the MPAP, 
whereby such trajectories are determined by the minimum of 
action, introduced as an effective cost function. In ifTsIl . such 
an approach has been implemented for an extension of the 
escape problem to an excitable FHN unit, having consid¬ 
ered the ghost separatrix as the terminating boundary. Due 
to fundamentally different terminating boundary conditions, 
one cannot apply the Hamiltonian theory per se for the prob¬ 
lem of first pulse emission, as we explain in more detail fur¬ 
ther below. Therefore, our approach cannot be interpreted in 
the genuine context of or be referred to as an extension of 
the Hamiltonian theory. In fact, the Hamiltonian system we 
consider should be interpreted as a set of effective equations 
integrated for relevant boundary conditions, whereby a partic¬ 
ular trajectory is singled out according to a certain predefined 
recipe. Our current goal is not to provide a systematic theory, 
but only to point out to striking similarity between the numeri¬ 
cally obtained MPAPs for the process of first pulse emission 
and the trajectories generated by the set of effective Hamil¬ 
tonian equations. This can be considered a preliminary stage 
of a study which may ultimately lead to derivation of an ap¬ 
propriate theory for the process of first pulse emission, whose 
basis would take into account certain aspects of the standard 
Hamiltonian formalism. 

In the remainder of this subsection, we briefly consider 
the main elements of the Hamiltonian approach, whose de¬ 
tailed description may be found in and clarify the dif¬ 
ferences emerging in case of the first pulse emission process. 
Within the Hamiltonian formalism, the optimal trajectories for 
a noise-driven escape process are obtained by variation-like 
approach which minimizes the appropriately defined ’’cost 
functional” S{J,t) along the set of possible activation tra¬ 
jectories J between the two fixed boundaries. In particular, 
the probability for noise to induce the Xi ^ Xf transition 
is given by p{xf\xi) = fj P[J]d{J}, such that d{N} de¬ 
notes integration along all the paths {57 = (xi,... ,X 7 v}} 
connecting Xi and x/. Each path is weighted by the probabil¬ 
ity P[J] oc where S is the cost function 

for the particular path. When D is small, the largest contribu¬ 
tion to p(xf\xi) comes from the path with the minimal cost 
function S'min = ram{Sj\J = (xi,... , xat}}. In other 
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words, small noise rarely gives rise to transition events, but 
once they occur, all the other orbits are suppressed in favor 
of the one corresponding to the minimal cost function. The 
transition probability then takes the asymptotic form 


p{xf\xi) = zexp[ = Xi.XN = Xf, (5) 


whereby the prefactor x is associated to ’’degeneracy” of the 
minimum, or rather the number of trajectories close to the one 
given by Smin I0- 

The form of S is selected so to reflect on one hand the im¬ 
pact of noise, whereas on the other hand to explicitly incor¬ 
porate the constraints between the variables and the stochastic 
terms which derive from the equations of motion. The nec¬ 
essary conditions for the minimum of such a constraint prob¬ 
lem can then be treated by the Lagrange multiplier technique. 
Ultimately, one arrives at a set of equations that may be in¬ 
terpreted as a Hamiltonian system, where the Lagrange multi¬ 
plier A enact the momenta conjugate to the system variables. 
In particular, for system O, the corresponding Hamiltonian 
set including the generalized momenta is given by 


dx o 

— = x-x /3 - y + 

— =e{X + b)+ TyPy 

dpx /, 2 \ 

— = -(1 - X )p^ - epy 

dpy 

dt 


= Px, 


(6) 


where p^ and py are the components of the momentum, while 
Xx and Ty represent the scaled noise intensities. By the lat¬ 
ter notation, Di and D 2 are conveniently expressed in terms 
of a single amplitude D (Di = XxD and D 2 = VyD). Sys¬ 
tem (01 has the Wentzel-Freidlin Hamiltonian H =Px{x — 
x^/3 — y) + Pye{x -f 6) -f ^pl + ^P^, whereas the tra¬ 
jectories connecting the initial state i and final state / are 
characterized by the action S = dt^{rxP^ + XyPy). At 
variance with system ([T]i, the fixed point of the system (01, 
{xo,yQ,Px,o,Py,o) = {-b,-b -f 6^/3, 0,0), is unstable for 
b > 1, which is corroborated by the positive real part of the 

characteristic exponents pi 2 = ^ ^^^ 2 ^ ^ ^ 

Conceptually, the Hamiltonian formalism is typically ap¬ 
plied to escape problems, where the stable equilibrium coex¬ 
ists with a certain saddle state, most often the saddle cycle 
dUl]. The aim is then to obtain the heteroclinic trajectory in 
the extended space, which emanates from the unstable man¬ 
ifold of the fixed point at f —?> —00 and reaches the stable 
manifold of the saddle cycle asymptotically at t ^ 00 , as 
well as tangentially (p —?> 0) for f 0. In principle, depend¬ 
ing on p, the dynamics in the extended space may also support 
the trajectories that do not settle at, but run across the cycle, 
as well as the trajectories that are repelled by the cycle, thus 
being reflected back to initial state. In ifTsIl . an extension of 
the escape problem to a single FHN unit has been consid¬ 
ered by utilizing the threshold-like behavior to construct the 
’’ghost” manifold (a separatrix in the asymptotic limit e —?> 0), 


which plays the role of a saddle structure. For such a sce¬ 
nario, the Hamiltonian approach yields the optimal trajectory 
connecting the unstable manifold of the fixed point to the sta¬ 
ble manifold of the ghost separatrix. 

Compared to this, the problem of first pulse emission is fun¬ 
damentally different because one should look for trajectories 
that cross the ghost separatrix, viz. leave the basin of attrac¬ 
tion of the fixed point. While the action surface in vicinity 
of the basin boundary should have one clearly defined global 
minimum, the physical picture near the spiking branch of the 
cc-nullcline typically involves multiple local minima of the ac¬ 
tion surface. A deeper theoretical analysis of these issues is 
beyond the scope of the current paper. The point we make 
here is as follows. Let us consider the trajectories that cor¬ 
respond to local minima of the action surface and adopt as a 
rule to select the solution that has the smallest momentum at 
the terminating boundary. Then, for different system configu¬ 
rations, the comparison of the MPAPs numerically obtained 
from an ensemble of stochastic realizations reveals significant 
similarity to the respective trajectories generated by the set of 
effective Hamiltonian equations. 

In terms of numerical treatment, obtaining the relevant tra¬ 
jectories from the effective system (0) comprises a bound¬ 
ary value problem. At the initial moment ti, the coordinates 
{x(ti), y{ti),px{ti),py{ti)) lying on the unstable manifold of 
the fixed point are specified according to the method provided 
in ll2ll] . which connects the coordinates in the configuration 
space with the components of the momentum. By this method, 
it follows that the initial coordinates in the extended space can 
be parameterized via a single angular variable (p G (0, 27r). 
The different trajectories are then obtained by sampling 1000 
equally spaced initial conditions that cover the entire range 
of (p values. The moment of arriving at the second bound¬ 
ary ts, as well as the coordinates {x{ts), y{ts),Px(ts),Py{ts)) 
are not explicitly specified. In effect, we apply a shooting ap¬ 
proach that involves a set of different trajectories with partic¬ 
ular initial conditions, whereby all the trajectories terminate 
at the spiking branch of the cubic nullcline. System (0i is in¬ 
tegrated by a standard solver implementing the fourth-order 
Runge-Kutta routine. The cost function S is calculated along 
each of the trajectories, and we single out the trajectory that 
corresponds to a local minimum having the smallest value of 
angular momentum at the terminating boundary (typically of 
the order of 10“'* or less). Note that the analogous numeri¬ 
cal approach is applied in case of two interacting units. The 
only difference is that the initial conditions are parameterized 
in terms of three, instead of a single angular variable. 


B. Examples of MPAPs and the method’s persistence under 
increasing noise 

We systematically examined how the MPAPs change un¬ 
der variation of the ratio of external vs internal noise intensity 
fxl'Ty It has been established that the MPAP profiles are 
not sensitive to gradual changes over the whole range of Vxlxy 
values. In particular, when slowly increasing rx/xy, the tra¬ 
jectories are found to either exhibit barely visible changes in 
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shape, or to converge to each other once the ’’ghost separatrix” 
is crossed. Therefore, in qualitative terms one may single out 
three characteristic forms of solutions, corresponding to cases 
{ra:,ry) = {l,0),{r^,ry) = (0,1) and {r^.Ty) = (1,1), 
which can be held representative for the problem of first pulse 
emission. In other words, the topological features of MPAPs 
are primarily sensitive to whether a particular noise source or 
both sources are present in the system O- 

The other point one has to consider is the impact of the 
noise intensity D. Naturally, the physical picture described 
above is maintained for sufficiently small D. The analysis 
provided in subsection IIIIDI will show that the term ’’suffi¬ 
ciently small” D effectively implies that the system lies suffi¬ 
ciently away from the stochastic bifurcation underlying transi¬ 
tion from stochastically stable fixed point to the noise-induced 
oscillations. Above the stochastic bifurcation, the attractive 
power of the fixed point is no longer felt, such that the noise 
can move the phase point away from equilibrium without an 
opposing force. Thus, if one focuses on MPAPs for any of 
the three characteristic Vx/ry ratios, the impact of noise will 
become substantial as one approaches the stochastic bifurca¬ 
tion, and will become overwhelming above the stochastic bi¬ 
furcation. Note that for large D, the MPAPs determined via 
the method introduced in subsection llll Al lose physical mean¬ 
ing, because the fluctuations over an ensemble of stochastic 
realizations grow too large to be accurately described by the 
maximums of the prehistory probability density H{x,y,t)- 
The influence of noise is reflected in the increase of the skew¬ 
ness and the kurtosis of the distribution of system’s variables 
for different stochastic realizations at arbitrary t. 

Let us now consider some examples of MPAPs in or¬ 
der to corroborate the points stated above. First we illustrate 
the validity of the method used to obtain the MPAPs, see 
Fig. 12 a). The figure refers to the case {rx,ry) = (0,1) for 
D = 0.00003, the noise value substantially below the stochas¬ 
tic bifurcation. By taking ten arbitrary realizations of the first 
pulse emission process, it is shown that the stochastic realiza¬ 
tions indeed cluster around the MPAP, indicated by the open 
circles. We have observed that the distribution of (x, y) values 
at fixed t for different stochastic realization is narrow around 
the terminating boundary, and broadens towards the initia¬ 
tion point. With increasing noise intensity D, the changes 
in the shape of the MPAPs become clearly visible around 
D Ri 0.00015, the value close to the onset of stochastic bifur¬ 
cation, cf. Fig. I2b). The effect of noise is felt particularly 
strong in the region of phase space before the ghost- separa¬ 
trix, which is indicated by the dotted line. An interesting point 
is that the larger noise may also have an inhibitory effect on 
the process of first pulse emission, in a sense that the phase 
point which has already crossed the ghost separatrix may still 
diffuse back to the basin of attraction of the fixed point. 

We further highlight how the topological properties of 
MPAPs depend on the characteristic ratio Txjry. As an¬ 
nounced earlier, the discussion is put into a broader per¬ 
spective by comparing the MPAPs to the trajectories gen¬ 
erated by the effective Hamiltonian system (|6]l according to 
the prescription provided in Subsection llll Al Taking the case 
Tx ■ ry = 0 : 1 as an example, we first demonstrate how the 



FIG. 2. (Color online) MPAPs for an excitable FHN unit, (a) 
The main frame shows that the activation paths for different stochas¬ 
tic realizations (solid lines) cluster around the MPAP (open blue 
circles). The ratio of external vs internal noise intensities is '■ 
Vy = 0 : 1, whereas the noise intensity is D = 0.00003. The 
a;-nullcline {NULL) and the canard-like trajectory (CNRD) per¬ 
taining the ghost separatrix are shown by the dashed and the dotted 
line, respectively. In the inset is displayed the enlarged view of the 
region of phase space before the CNRD. In (b) is illustrated how 
the MPAP profile changes under increasing D for the fixed ratio 
rx ■ Ty = 0 ■. 1. The MPAPs shown correspond to D = 0.00003 
(circles), D = 0.00015 (squares) and D = 0.001 (diamonds). These 
noise intensities lie substantially below, in vicinity and above the 
stochastic bifurcation, respectively. 


problem of first pulse emission is different from the extension 
of the escape problem to an excitable FHN unit addressed in 
Q. Apart from the ghost separatrix, denoted by the dotted 
line. Fig. [2a) also shows the optimal trajectory for the escape 
problem obtained using the standard Hamiltonian formalism. 
The latter trajectory is indicated by the solid line which ap¬ 
proaches, but does not cross the ghost separatrix, because the 
optimal escape path cannot intersect the system’s manifold. 
The MPAP numerically determined for the process of first 
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----NULL 

CNRD 
• (D,, D^)=(0.00009, 0) 
■ D=D=0.00009 


FIG. 3. (Color online) Extended analysis of MPAPs. (a) is in¬ 
tended to illustrate the difference between the first pulse emission 
problem and the generalized escape problem for a FHN unit. The 
setup and the style of presentation are the same as in Fig. |3a). but 
two additional curves are presented. The solid black line approach¬ 
ing the CNRD indicates the optimal trajectory obtained for the es¬ 
cape problem via the standard Hamiltonian approach, whereas the 
dash-dotted line denotes the trajectory generated by the system © 
according to the recipe described in subsection IIII Al (b) shows the 
MPAPs for '■ Ty = 1 ■. 0 (circles) and ■ Vy = 1 ■. 1 (squares), 
together with the corresponding trajectories (dash-dotted lines) gen¬ 
erated by the system ®. In both instances, the noise intensity is 
D = 0.00009. 


pulse emission is represented by the open squares. The par¬ 
ticular MPAP is obtained for small, but finite noise D. Such 
D values may actually be referred to as intermediate, because 
on one hand, they lie below the critical domain giving rise to 
stochastic bifurcation, but on the other hand, they cannot be 
considered as asymptotically small noise [D 0) where the 
theory of large fluctuations would naturally apply. It can be 
seen that the initial part of the MPAP matches very closely 


the optimal trajectory calculated by the standard Hamiltonian 
formalism, but also departs from it well before reaching the 
ghost separatrix. Nevertheless, an interesting finding is that 
the trajectory generated from the effective system (|6]l, cf. the 
thick solid line that crosses the ghost separatrix, matches quite 
closely the MPAP along the entire trajectory relevant for 
the process of first pulse emission. This point suggests that 
a theory using certain aspects of the standard Hamiltonian ap¬ 
proach could potentially be derived to characterize the first 
pulse emission process. 

Fig. |3!b) is intended to compare the MPAPs for the two 
remaining characteristic Tx/vy ratios at D values substantially 
below the stochastic bifurcation. As expected, the respective 
trajectories show significant differences well before crossing 
the ghost separatrix. The trajectories generated by the effec¬ 
tive Hamiltonian equations in the fashion described in Sub¬ 
section |IIl3 again seem to closely match the MPAPs ob¬ 
tained by averaging over the ensemble of stochastic realiza¬ 
tions. The profile of the escape paths in vicinity of equilib¬ 
rium in Fig. [2b) suggests that the noise-induced linearization 
ll^ [^ takes place in presence of the external noise. In par¬ 
ticular, the latter is found to smear the effect of nonlinearity, 
such that the threshold-like behavior becomes smoother. 


C. An insight into the activation processes driven hy external 
or internal noise 

Before proceeding to numerical results regarding the statis¬ 
tical properties of the activation process, let us briefly consider 
the conceptual differences between the cases where noise in¬ 
fluences the dynamics on the fast {Di > 0, D 2 = 0) or the 
slow characteristic timescale {Di = 0,D2 > 0). Note that 
the ’’mean activation times” obtained from approximations in¬ 
troduced here are not intended to be compared quantitatively 
with the actual stochastic averages from Sec. IIIIDI but are 
rather aimed at gaining qualitative insight into the distinct 
mechanisms by which the two noise sources affect the activa¬ 
tion process. To this end, we use the standard approach which 
consists in reducing the original dynamics, given by dTJ with 
Di or D 2 set to zero, to an appropriate Langevin equation of 
the form dz = —U'{z)dt -F y/2DdW(t), where U(z) denotes 
the effective potential. 

If only Di is present, one may conveniently exploit the 
sharp separation between the two characteristic time scales. 
The analysis is confined solely to the fast variable subsystem, 
which is first rewritten as 

dx = - + y/^idW, (7) 

ox 

with U{x,y) = —^x"^ + ■^x'^ + xy. In the last expression 
for U, y may simply be seen as a parameter, because the 
evolution of the y-variable takes place on a timescale much 
slower than that of the x-variable ll44l] . Another useful point 
is that U has the form of a double-well potential. Its two lo¬ 
cal minima, as well as the local maximum, correspond to x- 
values which are the roots of the equation for the x-nullcline 
X — ^x^ — y = 0. In particular, the local minima are given by 
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FIG. 4. (Color online) Assessing the impact of two distinct noise 
sources on the activation process of a single FHN unit. The curves 
refer to approximate dependencies of ’’activation times” T on Di 
(dashed line) and D 2 (solid line). Note that the introduced approx¬ 
imations are crude, so that the results can only be considered in a 
qualitative fashion. Still, the two curves accurately predict that the 
activation process led by D 2 is comparably faster than the one led by 
T>i. 


the solutions xi and X 3 lying on the refractory and the spik¬ 
ing branch, respectively, while the local maximum coincides 
with the X 2 solution on the unstable branch of the x-nullcline 
(xi(y) < X 2 (y) < xsly)). Then, the activation process may 
be interpreted as a jump from the well at the refractory branch 
to the one at the spiking branch, whereby the phase point has 
to overcome the potential barrier provided by the local maxi¬ 
mum. 

Given that U has a double-well shape, the escape time for 
a particle to jump from xi(y) to a point near X 3 (y) can be 
approximated by the well-known Kramers formula jd^l - ld^] 
T = case, one 

y/\U"{x2)\U’’{xi) 

may fix the y = —2/3 value, which corresponds to the min¬ 
imum of the cubic nullcline. Having determined the appro¬ 
priate a;i(—1) and X 2 {— 1 ), we have obtained the dependence 
T{Di) shown by the dashed line in Fig. |4l Naturally, given 
the crudeness of the introduced approximations, both in terms 
of system dynamics and the boundary conditions, the obtained 
T values should not be compared to the mean activation times 
for our activation problem. 

The scenario where only D 2 is present cannot be ap¬ 
proached the same way as above, because y can no longer be 
treated as a parameter. Still, one may implement the adiabatic 
elimination of the fast variable to find the effective potential 
that influences the y dynamics. Though the roots of the equa¬ 
tion describing the slow manifold y = x — ^x^ may be used 
explicitly, the drawback is that such solutions are not easily 
handled analytically. In order to keep the subsequent form 
of effective potential analytically tractable, it is more conve¬ 
nient to approximate the relation between y and x in vicinity 
of the minimum of the cubic nullcline by a simple relation 


y{.x) = ym + - XmY, where {xm,ym) = (-1, -2/3) 

refer to coordinates of the minimum. In other words, in prox¬ 
imity of the minimum, the cubic dependence has been re¬ 
placed by a quadratic approximation. Note that k — —2xm 
is determined by taking the second derivative of the equation 
for the x-nullcline. From the expression for y{x), one readily 
finds that x = Xm ^ y/y — ym, whereby the plus sign solu¬ 
tion is relevant for our activation problem. It follows that the 
equation for the dynamics of y may be written as 

dy = e{xm + \/y - ym + h)dt + \/2D^dW, (8) 

such that the corresponding effective potential reads U = 
e{xm + b)y + |e(y - ym)^/^. 

Given that U is not a double-well potential, one cannot use 
the Kramers-like equation for the mean first-exit time. In¬ 
stead, it is appropriate to use the general form 

-| /*a pu 

T=4- (9) 

^2 Ji Jr 

derived from the Fokker-Planck approximation to a problem 
involving a single absorbing boundary a and a single reflect¬ 
ing boundary r Note that i in one of the integration 

limits refers to the injection point (initial location), which in 
our activation problem corresponds to the deterministic fixed 
point {xeq, yeq) = (—&, —b + b^/3). Regarding the reflecting 
boundary, it has often been found that the final results are not 
affected by its particular value 1110, so r may readily be 
set to r ——00 to simplify the calculations. As for the ab¬ 
sorbing boundary, it generally concerns the terminating point 
of the activation path and remains a free parameter that can¬ 
not be known a priori. Nevertheless, we may use the results 
for the MPAP solutions from Sec. IIIIBI and simply read the 
necessary coordinates of the terminating point (y = —0.308), 
which illustrates the complementary nature of the methods 
applied. Note that the integrals such as the ones in (|9|l are 
typically resolved by introducing convenient approximations. 
In the particular case, we have used the approximate solution 
provided in fill] . The obtained curve T{D 2 ) is shown by the 
solid line in Fig. |4] 

It has already explained that the results in Fig. |4] cannot 
be considered reliable in quantitative sense given the crude¬ 
ness of the approximations involved. Still, certain qualitative 
insight have been gained. For instance, the analysis above in¬ 
dicates that the activation processes led by Di or D 2 have two 
considerably distinct backgrounds, whereby only the mecha¬ 
nism of the former may be interpreted in analogy to a jump 
over the barrier that separates two potential wells. Further, the 
curves in Fig. |4]turn out to be consistent with the general trend 
for a single FHN unit demonstrated in the next subsection, 
according to which the activation is more easily excited by D 2 
than Di. 

D. Statistical features of the activation process 

The statistics of activation events is characterized in terms 
of the mean TFP t{Di,D 2 ) and the associated coeffi¬ 
cient of variation R{Di, 02 ). The former is an average 
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FIG. 5. (Color online) Statistical features of activation process influ¬ 
enced by Di and D 2 - In (a) is displayed how the mean TFPs r aver¬ 
aged over an ensemble of different stochastic realizations depend on 
D\ and D 2 , whereas (b) concerns the associated coefficient of varia¬ 
tion R{Di, D 2 ). 'r{Di, D 2 ) exhibits three characteristic regimes of 
behavior, whereby the transition from the large values to the plateau 
region is found to qualitatively correspond to stochastic bifurcation 
from the stochastically stable fixed point to continuous oscillations. 
An indication on the noise intensities that give rise to stochastic bi¬ 
furcation is provided in the inset of (a). The latter shows bifurcation 
curve D 2 {Di) obtained for the approximate model dl lb of the origi¬ 
nal stochastic system 0. 


of TFPs for different stochastic realizations t{Di,D 2 ) = 

Ur 

— ^ Ti{Di, 02 ). The stochastic paths taken into account 

i—l 

satisfy the specified boundary conditions, such that the trajec¬ 
tories emanate from the deterministic fixed point and termi¬ 
nate at the spiking branch of the cubic nullcline. The coef¬ 
ficient of variation is defined as the normalized variation of 
activation times 


R{Di,D2) = ( 10 ) 

{T-i) 

where (•) refers to averaging over an ensemble of stochastic 
realizations. Quantity R is intended to describe the regular¬ 


ity of the activation process, in a sense that the smaller R 
indicates that the TFPs deviate less from the mean value. 
Note that numerical simulations for all the considered setups 
are carried out by implementing the Euler integration scheme 
with the fixed time step 6t = 0.002, having verified that no 
changes occur for smaller 6t. All the results for the mean 
TFPs and the associated variances are obtained by averaging 
over an ensemble of 5000 different stochastic realizations of 
the activation process. In each realization, the initial condi¬ 
tions of a unit coincide with the deterministic fixed point. 

The fields t{D 1 , 02 ) and R{Di, D 2 ) for a single excitable 
unit are plotted in Fig. |3a) and Fig. lUb), respectively. The 
obtained profiles indicate that the mean TFPs are more sen¬ 
sitive to variation of Di, whereas R shows strong dependence 
on Z? 2 - Regarding the mean TFPs, one may distinguish three 
characteristic regimes, including (i) long TFPs, encountered 
for small Di and D 2 , (ii) the plateau region, comprising in¬ 
termediate Di and intermediate to large D 2 values, as well as 
(in) short TFPs, found for large Di irrespective of £> 2 . 

Given that the considered stochastic process is influenced 
by two sources of noise, one finds that the transitions between 
the different regimes are gradual rather than sharp, whereby 
the ’’boundaries” are naturally smeared due to action of noise. 
Nevertheless, in terms of theory, it is reasonable and well jus¬ 
tified to discuss the physical background giving rise to transi¬ 
tions between the different regimes. What we postulate is that 
the transition between the domains (i) and (ii) may qualita¬ 
tively be accounted for by the fact that the excitable unit un¬ 
dergoes stochastic bifurcation induced by Di and D 2 . Note 
that the phenomenological stochastic bifurcation we 

refer to corresponds to the noise-induced transition from the 
stochastically stable fixed point (stationary probability distri¬ 
bution P{x, y) focused around the fixed point) to the stochas¬ 
tically stable limit cycle (stationary probability distribution 
P{x,y) showing non-negligible contribution for {x,y) val¬ 
ues along the spiking and the refractory branches of the x- 
nullcline). Intuitively, one understands that the fixed point can 
be considered stochastically stable if the amplitude of fluctu¬ 
ations around the fixed point is of the order of noise intensity. 
By the same token, one may perceive a limit cycle as stochas¬ 
tically stable if the general structure involving two branches 
of the x-nullcline (the spiking and the refractory branch) is 
preserved under the action of noise. 

In support of associating the transition between the domains 
of large TFPs and the plateau region with the stochastic bi¬ 
furcation, one may invoke the qualitative argument that above 
the bifurcation, the attractive power of the fixed point is ef¬ 
fectively no longer felt, in a sense that noise can drive the 
phase point away from equilibrium without meeting a strong 
opposing force. At the level of mean TFPs, this should be re¬ 
flected as follows. Below the stochastic bifurcation, the mean 
TFPs are expected to be longer, while above the bifurcation 
they should substantially reduce and also become fairly in¬ 
sensitive to further increase of noise. In other words, once the 
fixed point becomes stochastically unstable, the gross effect 
of noise is the same because the fixed point holds no attractive 
power to resist its action. 

Having established that the boundary between the domain 
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of long TFPs and the plateau region in Fig. |5la) should co¬ 
incide with the values that give rise to stochastic 

bifurcation, our next goal is to try to obtain these values an¬ 
alytically. To do so, we derive a deterministic model based 
on Gaussian approximation of the original system ([T]l- Ac¬ 
cording to Gaussian approximation, all the cumulants above 
the second order are assumed to vanish ll36l - [^ . One is ul¬ 
timately left with a set of five equations describing the dy¬ 
namics of the first moments mx{t) = E[x{t)] and -iriy = 
E[y{t)], the variances Sx{t) = E[{x{t) — mx{t)y] and 
s„(t) = E\(y(t) — as well as the covariance u(t) = 

E[{x{t) - mx{t)){y{t) - my{t))]: 

1 3 

mx = nrix- -m^ - rrixSx - rriy 
niy = e{mx + b) 

s'x = 2sx{l — rri^ — Sx) — 2u + 2Di 
Sy = 2 cu -\- 2D2 

ii = u{l - ml - Sx) +esx - Sy. ( 11 ) 

Note that the impact of noise in (fTTT) is described only by the 
respective intensities Di and D 2 , which can then be treated as 
bifurcation parameters. In particular, the bifurcation analysis 
shows that the approximate model (fTTl i displays direct super¬ 
critical Hopf bifurcation, which qualitatively accounts for the 
stochastic bifurcation exhibited by the excitable unit ([T]i- The 
obtained bifurcation curve D 2 {Di), plotted in the inset of Fig. 
|3a), and hence the transition from the large values of mean 
TFPs to the plateau region. 

Better understanding of the nature of the activation process 
behind the r and R dependencies from Fig. |3a) and Fig. |5fb) 
may be gained by examining how the distribution of TFPs 
over an ensemble of different stochastic realizations changes 
under variation of Di and 02- Let us note first that even for 
a typical escape problem, such an issue has rarely been ad¬ 
dressed in the literature and is difficult to approach analyti¬ 
cally. In particular, for the distribution of first exit times in 
a typical escape problem, the only rigorous result so far has 
been found within the framework of large fluctuations theory 
S. It states that for the sufficiently small noise intensity, 
the distribution of the first exit times asymptotically acquires 
exponential form. If interpolated to the case of an excitable 
FHN unit, one should expect the latter statement to apply 
regardless of whether noise is added solely to the fast or the 
slow variable. Still, note that the mentioned result could con¬ 
cern only the extension of the escape problem to excitable sys¬ 
tems, where the terminating boundary set would be given by 
the ’’ghost separatrix”. Recall that we have already explained 
the conceptual difference between such a problem and the ac¬ 
tivation problem we consider, where the focus exclusively lies 
with the spiking response of an excitable unit. 

The characteristic examples illustrated in Fig. |6] suggest 
that the activation process dominated by D 2 gives rise to the 
exponential distribution of TFPs, whereas the distributions 
generated by the prevailing Di conform to Lorentzian-like 
profile with the cut-off at small TFP values. Regarding the 
former, one should emphasize that the exponential distribu¬ 
tion of the inter-event intervals is typically associated with 



FIG. 6. (Color online) Impact of Di and D 2 on the distribution of 
TFPs P{t) obtained for an ensemble of different stochastic real¬ 
izations. D 2 alone typically gives rise to an exponential distribu¬ 
tion of TFPs, which is consistent with the Poisson process. Un¬ 
der prevailing D\, the distributions are found to conform to a uni- 
modal, Lorentzian-like profile. The displayed examples refer to cases 
Di = 0.0007, D 2 = 0 (solid squares), Di — 0.0001, D 2 ~ 0.0001 
(solid triangles), Di = 0.02, D 2 = 0 (solid circles), Di — 0, D 2 = 
0.0001 (empty circles) and Di = 0, D 2 = 0.02 (empty squares). 


the Poisson process iH, and it is indeed not unexpected to 
find an excitable unit under small amount of noise to act as a 
Poisson generator lEl] . Nevertheless, an interesting finding is 
that the activation events led by Di seem to be derived from 
some process other than the Poissonian, though the asymp¬ 
totic distribution of TFPs is still exponential. Note that the 
qualitative distinction between the average effects of Di and 
D 2 has already been commented on in subsection III Al The 
remark on the two distribution types holds if the leading noise 
term is much stronger than the remaining one. Nevertheless, 
while the validity of this statement is maintained even under 
substantial increase of Di as the leading term, the analogous 
point for D 2 applies up to a certain value, as explained below. 

For {Di,D 2 ) corresponding to the plateau from Fig. 
13a), one generally encounters the exponential distribution of 
TFPs or some of its modifications. In this context, note that 
the increase of R in Fig. |3b) remains fairly slow, if any, for 
D 2 values that warrant R < 1, but becomes steep once R 
passes 1 around D 2 ^ 10“^. The approximate boundary at 
i? = 1 coincides with the coefficient of variation one obtains 
for the exponential distribution of TFPs. It is further found 
that the sufficiently large D 2 values where i? > 1 give rise to a 
peculiar regime where the ensemble of activation events splits 
in two sharply distinct classes, such that the one with the small 
TFPs dominates, but the contribution from the events with 
large TFPs is non-negligible. The examples of stochastic ac¬ 
tivation paths for large D 2 , including those with long TFPs, 
where the phase point rebounds onto the refractory branch of 
the a;-nullcline before the transition to spiking branch is trig¬ 
gered, are already provided in Fig. [3b). Note that the inten- 












tion there has been to illustrate the setup for which the analyt¬ 
ical method of calculating the MPAPs fails. 


12 


IV. CASE OF TWO UNITS: MPAPS AND STATISTICAL 
PROPERTIES OF THE ACTIVATION PROCESS 


In this Section, we investigate how the form of coupling af¬ 
fects the first pulse emission process for two units subjected 
to external and internal noise. The first subsection is focused 
on the MPAPs obtained for the setups with linear or nonlin¬ 
ear interactions, whereas the second subsection concerns the 
statistical features of the activation process. 


A. Dependence of MPAPs on the form of coupling 

The MPAPs are obtained by the numerical method pre¬ 
sented in subsection IIII Al Nevertheless, before proceeding 
with the results, we first consider what counts as a two- 
unit activation event by providing the appropriate boundary 
conditions. In terms of the initial conditions, the activation 
paths of units described by (l2]i or Q begin at the equilibrium 
(^l,eg ; yi,eg; ^2,eg ; y2,eg) ( bb /3, 6, bb /3). 
The terminating boundary conditions are specified in such a 
way that the phase points of both units should reach the spik¬ 
ing branch of their respective Xi-nullclines. This definition 
implies that the time-to-first pulse for a couple of units is de¬ 
termined by the slower-firing unit. Note that the profile of the 
spiking branch for a coupled unit is quite similar to that of a 
single unit, which can readily be verified. In particular, for the 
given unit one may approximate the a:-coordinate of the other 
unit within the coupling term by its value at the deterministic 
fixed point, around which the actual values fluctuate for most 
of the time. 

Regarding the coupling strengths, note that we are inter¬ 
ested only in c values that are subcritical with respect to Hopf 
bifurcation, viz. c is always set so that the deterministic ver¬ 
sions of dU or (O admit the locally stable fixed point. This is 
consistent with the goal to examine the noise-driven first pulse 
emission process and the fashion in which it is modified by the 
presence of coupling. The fact that the c values are subcritical 
also implies that in the deterministic case, the large excitation 
of one unit (initial conditions far from equilibrium) cannot in¬ 
duce a spiking response of the other unit which lies at equi¬ 
librium. One should point out that the adopted formulation 
of the activation problem for two units is by far more appro¬ 
priate than the alternatives, including the one cast in terms of 
the average variables (xi X2)/2 and {yi y2)/2. Without 
stating the details, we make a remark that the latter approach 
would have several drawbacks. On one hand, the dynamics 
associated to synchronization would smear the physical pic¬ 
ture relevant for the activation problem, whereas on the other 
hand, there would be no immediate generalization from the 
case of a single unit to the setups with two units. 

Now let us focus on the scenario where two excitable units 
interact via linear couplings. The symmetry of interactions is 
reflected in the point that the MPAP trajectories are identical 
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FIG. 7. (Color online) MPAPs for two coupled units. In (a) are 
displayed the MPAPs for two units (open and solid squares) in¬ 
teracting via linear couplings of subcritical strength c = 0.04. The 
solid lines denote ten sample paths corresponding to the process of 
first pulse emission. The results are obtained for ra, : r,, = 1 : 1 
and D = 0.00005. The dash-dotted line indicates the trajectory 
generated by the extended system GIJ. In the inset are shown 
the MPAPs for the same characteristic ratio and noise intensities 
D = 0.00002 (squares) and D = 0.0005 (circles), (b) provides 
a comparison between the MPAPs obtained for linear (open and 
solid circles) and nonlinear interactions (open and solid diamonds). 
In the latter case, the subcritical coupling strength is c = 0.07, while 
the noise parameters are the same as in (a). In the inset is shown the 
enlarged view of the portion of phase space before the CNRD. 


for both units. This is corroborated in Fig. Ha), which shows 
as an example the respective MPAPs of the two units (solid 
and open squares) for the characteristic ratio : ry = 1 : 1 
at small D substantially below the stochastic bifurcation. Sev¬ 
eral realizations of the first pulse emission process are plotted 
to demonstrate the clustering around the MPAP. The ap¬ 
proximate matching between the units’ most likely trajecto- 
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lies is found to persist for all {rx,ry) characteristic setups. 
The changes of the MPAP profile under increasing noise in¬ 
tensity D are illustrated in the inset of Fig. EJa). It has already 
been mentioned that depending on c, the units coupled in a 
linear fashion may display two different regimes where the 
fixed point is stable. These regimes have been referred to as 
excitability proper and generalized excitability, whereby the 
latter involves coexistence between the stable fixed point and 
the stable limit cycle created in a global fold-cycle bifurcation. 
Comparing the cases where c is sub- or supercritical with re¬ 
spect to global bifurcation, we have established that for small 
but finite D, the profile of MPAPs for all three characteristic 
ratios Tx/vy does not appear to show significant differences. 

As in case of a single unit, we make a brief remark regard¬ 
ing the trajectories generated by the effective set of Hamilto¬ 
nian equations. Compared to (|6]l, the equations for the dynam¬ 
ics of the extended system are modified to include the interac¬ 
tion terms 


dxi 

dt 


= Xi- x\li 


Vi + TxVx.i + c{x^ 


dVi 

dt 

dpx,i 

dt 

dpy^i 

dt 


= e{xi -I- 6) -I- XyPy^i 
— (1 )Px,i ^Py,i ^Px,j 

Px,i t 


Xj) 


( 12 ) 


whereby i,j G {l,2},i ^ j denote the units’ indices. The 
numerical treatment of the above system again involves in¬ 
tegration for the initial conditions set on the unstable mani¬ 
fold of the saddle point. In particular, in configuration sub¬ 
space one chooses the initial values {xi, pi) that lie on a four¬ 
dimensional sphere of a very small radius, which encloses 
the deterministic fixed point J/i.eq; a; 2 ,e(j: J/ 2 ,eq)- Nat¬ 

urally, the points located on a four dimensional sphere are 
parametrized with three independent angular variables. The 
initial values of the generalized momenta are then obtained 
using the prescription provided in ll^ . An interesting finding 
is that the trajectories of the extended system, selected by the 
rule described in subsection llll Al again closely match the nu¬ 
merically obtained MPAPs, cf. the dotted line in Fig. Ha). 
This striking similarity further evinces that the theory in the 
spirit of the Hamiltonian approach may potentially be derived 
for the problem of first pulse emission. 

An issue that requires to be addressed is whether and how 
sensitive are the topological features of the MPAPs with re¬ 
spect to linear/nonlinear form of coupling. The comparison is 
facilitated in Fig. |2lb), where the corresponding MPAPs for 
the scenarios with linear (circles) and nonlinear interactions 
(diamonds) are plotted together. The data shown are obtained 
for the same Di and D 2 , whereas the coupling strengths are 
analogous in terms of distance from the Hopf bifurcation. The 
trajectories corresponding to the different units for the same 
system configuration are distinguished by the solid and open 
symbols. Note that in case of nonlinear interactions, there 
is no symmetry to warrant that the respective MPAPs of 
the two units would be identical for arbitrary values of ex¬ 
ternal and internal noise. It turns out that the MPAPs are 


approximately identical for intermediate D sufficiently below 
the stochastic bifurcation. Nevertheless, we have also found 
that the spread of the MPAPs in case of nonlinear couplings 
increases with D. As for the effects of linear vs nonlinear in¬ 
teractions, one observes significant differences in the profiles 
of MPAPs within the initiation region, which is shown en¬ 
larged in the inset of Fig. |7jb). Also note that the MPAPs of 
the interacting units substantially depart from what is found in 
case of a single unit, cf. Fig. |2]and Fig. [3] 


B. Statistical properties of the activation process for coupled 

units 

Here we consider two types of numerical results, one refer¬ 
ring to the two-unit activation events, and the other concern¬ 
ing the correlation between the activation events on individual 
units. Recall that the ’’compound” activation event for a cou¬ 
ple of units requires that the phase points of both units have 
reached the spiking branch of the appropriate cci-nullcline, 
consistent with the definition adopted in Sec. IIV Al Adher¬ 
ing to this, we have calculated the mean TFPs Tc{Di,D 2 ) 
and the associated coefficients of variation Rc{Di, D 2 ) for 
the scenarios with the linear or the nonlinear couplings, fur¬ 
ther examining how the results are changed under variation of 
c. Though c is always selected to lie below the Hopf thresh¬ 
old, in the latter context one may still distinguish between the 
strongly and the weakly subcritical regimes. 

In terms of whether the mere form of coupling affects the 
statistical features of the activation process, it may be shown 
that the fields Tc and Rc exhibit qualitatively analogous de¬ 
pendencies for the linear and the nonlinear interactions. It is 
further found that the behavior of mean TFPs is in several as¬ 
pects different to that of a single unit, whereas the correspond¬ 
ing coefficient of variation Rc is only marginally dependent 
on c, displaying the ’’universal” behavior qualitatively similar 
to the one in Fig. |5jb). As an example on how the prop¬ 
erties of activation process change with c, in Fig. |8ja) and 
Fig. Ob) are illustrated the respective fields Tc{Di,D 2 ) for 
the strongly and weakly subcritical regimes in case of the lin¬ 
ear coupling. Both plots corroborate the existence of the three 
typical regimes already indicated in reference to Fig. O^)^ 
though one no longer speaks of stochastic bifurcation in vicin¬ 
ity of the Hopf bifurcation controlled by b, but rather of the 
one induced by the coupling strength. 

Nevertheless, several differences due to presence of cou¬ 
pling should be noted. First, as c is increased, the mean TFPs 
for small noise intensities (below the stochastic bifurcation) 
become shorter, see Fig. |8jc). Also, at small noise intensi¬ 
ties the activation times Tc of the couple are reduced when 
compared to the case of a single unit. However, for intermedi¬ 
ate Di and D 2 corresponding to the plateau region, the mean 
TFPs for the pair are larger than those found for a single 
unit. More importantly, the average activation times in this 
region seem to increase as c rises, see Fig. Etd). One may in 
fact discern a general trend that the differences between the 
three characteristic regions are gradually washed out as c is 
enhanced, which is manifested even more once c is supercrit- 
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FIG. 8. (Color online) Statistical properties of activation process for two units interacting via linear couplings, (a) and (b) refer to mean 
TFPs Tc{Di, D2) in a strongly subcritical (c = 0.01) and a weakly subcritical regime (c = 0.04), respectively. Note that the definition of 
the two-unit activation event requires that the phase points of both units have reached the appropriate terminating boundary set. (c) and (d) 
show how Tc (main frames) and Rc (insets) behave under increasing c. The noise intensities (Di, D 2 ) = (0.00014, 0.0008) fixed in (c) are 
representative for the domain of large TFPs from (a), whereas (d) is obtained for (Di, D2) = (0.154, 0.0002), the values corresponding to 
the plateau region in (a). In (e) is demonstrated how the distributions of TFPs over different stochastic realizations vary with Di and D2- Note 
that the unimodal distribution profile is preferred over the exponential form. The data are obtained for Di = 0, D 2 — 0.0001 (solid triangles), 
Di — 0.0005, D 2 = 0 (solid squares), Di = 0.01, D 2 = 0 (solid circles), Di — 0, D 2 = 0.02 (empty circles) and Di = 0.005, D 2 = 0 
(empty squares). 


ical. This is to be expected in the latter case, given that the 
activation process then becomes primarily deterministic, viz. 
it is only perturbed by the noise terms. Note that in all the 
three characteristic regimes the corresponding coefficients of 
variation are seen to reduce with c, cf. the insets in Fig. [8])c) 
and Fig. Od). The final remark on the impact of coupling is 
that the region with small mean TFPs apparently decreases 
with c. 

We have further examined how the distribution of activa¬ 
tion events over different stochastic realizations depends on 
the form and strength of coupling. One may state that the 
general conclusions reached in case of a single unit, cf. Fig. 
|6] persist for the activation events of coupled units, though 
the physical picture is more perturbed for the scenario with 
the linear than the nonlinear interactions. Recall that for a 
single unit, only the external noise has been found to induce 
deviations from the typical exponential distribution of events. 
For the scenario involving the nonlinear coupling, there may 
be more noise domains admitting some distribution other than 
exponential, but the latter still constitutes the prevailing form 
of behavior. However, under linear interactions, it turns out 
that both Di and D 2 may give rise to a form of TFP distri¬ 


bution completely absent in case of a single unit. The partic¬ 
ular form is unimodal, and is numerically best approximated 
by the lognormal profile, cf. Fig. [8l)e). One cannot suspect 
on the type of activation process producing such a distribu¬ 
tion, though the explanation on why it is different from both 
the scenarios with a single unit and two units interacting in a 
nonlinear fashion likely has to take into account the existence 
of global bifurcation 1^ . 

As far as the effects of coupling strength are concerned, the 
TFP distributions typically become narrower as c approaches 
the critical value, viz. the tail at longer activation times is re¬ 
duced compared to an uncoupled unit for the same {Di, 02 ). 
This point holds independent on the linear/nonlinear form of 
coupling. Thus, one may state that the impact of stronger cou¬ 
pling is expectedly reflected in the decrease of fraction of the 
longer individual activation events. 

The final point we address concerns the relation between 
the individual activation processes on the units making up the 
pair. This issue may be approached from two angles, either 
by examining the mean difference in the single unit TFPs, 
or by analyzing the correlation between the individual acti¬ 
vation events. On the former point, we introduce a measure 
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FIG. 9. (Color online) Relationship between TFPs of individual units making up the pair. In (a) and (b) are plotted the differences in mean 
TFPs At{Di, D 2 ) and the associated coefficients of variation Rat{Di, D 2 ), respectively. The data are obtained for the linear couplings 
of strength c = 0.04, and similar results are found for the setup with nonlinear interactions, (c) and (d) show the correlation coefficients 
p{Di, D 2 ) between the TFPs for the different stochastic realizations in cases of linear (c = 0.04) and nonlinear interactions (c = 0.06), 
respectively. 


of coherence between the individual activation events aver¬ 
aged over an ensemble of rir different stochastic realizations, 

n-r 

At = — X] |Tfe 1 — Tfe_ 2 |j as well as the associated coeffi- 

k=l 

cient of variation Rat- At(Z9i, D 2 ) is intended to describe 
the net effect of how much the interaction is able to enforce 
matching between the activation times of noise-driven units. 
To understand the relevance of this, one should recall that the 
coherence of the first units’ responses is an issue quite distinct 
from that of asymptotic synchronization between the spiking 
series. 

In Fig. |9l)a) and Fig. |9jb) are plotted the dependencies 
At{Di, D 2 ) and Rat{Di, D 2 ) for the case of two units in¬ 
teracting via linear couplings, but qualitatively similar results 
are found for the nonlinear interactions as well. From Fig. 
|9l)a) one learns of a general tendency for the mean difference 
to increase with the internal noise. For smaller fixed D 2 , the 
spread of individual TFPs reduces under the action of exter¬ 
nal noise, but such an effect is lost once D 2 overwhelms the 
system dynamics. Note that for large internal noise At be¬ 
comes comparable to Tc for the coupled units, which suggests 
a typical scenario where the activation process of one unit is 
rapid, whereas the pulse emission of the other unit is substan¬ 
tially delayed. Naturally, for the {Di,D 2 ) domain support¬ 
ing small Tc, the At dependence tells us that both the units 


emit pulses virtually at the same time. Regarding Fig. |9jb), 
the gross effect is that the variation of the difference of single 
unit TFPs shows a steep increase above some ’’threshold” ex¬ 
ternal noise, whose value becomes larger as D 2 is enhanced. 
Similar dependence is found for the nonlinear interactions, but 
with the reverse roles of Di and 02 - 


Next we examine the correlation of individual activation 
times for units comprising the pair. The correlation is quanti¬ 
fied by the correlation coefficient between the times-to-first- 


pulses of single units p = 


{rk,irk,2)-{rk,l){'rk,2) 


v/V ('^k,2)-Pk,2P ’ 


where the angled brackets denote averaging over an ensem¬ 
ble of different stochastic realizations. The fields p{Di, D 2 ) 
for the setups with the linear and nonlinear interactions are 
shown in Fig. |9jc) and Fig. |9jd), respectively. A common 
ingredient in both cases is that there exists a certain value of 
D 2 , above which the noise effects prevail. There the first-time 
pulses are neither correlated nor anti-correlated, viz. the cor¬ 
relation coefficient lies around zero. However, for smaller D 2 , 
the correlation between the individual activation events sub¬ 
stantially depends on the form of coupling. For linear interac¬ 
tions, small Di then facilitates strong correlation, whereas for 
larger Di the TFPs become significantly anti-correlated. On 
the other hand, for the nonlinear couplings p displays a more 
homogeneous dependence on Di and 02- In fact, the corre¬ 
lation is strong for small noise intensities, and it reduces with 
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the increase of both Di and D 2 , whereby the rate of decline 
depends more on D 2 than Di. 

As an interesting remark for future study, we indicate a po¬ 
tential link between the form of correlation of the single unit 
TFPs and the synchronization state of the stochastic units’ 
time series ll45[ - i^ under the given parameter set. In par¬ 
ticular, note than in case of linear interactions the time se¬ 
ries of two units display a constant phase shift for small Di, 
which coincides with the domain where the T FPs are corre¬ 
lated. Nonetheless, the time series show phase slips and am¬ 
plitude fluctuations for larger D^, where the activation times 
are uncotTelated or anti-correlated. In a similar fashion, for the 
nonlinear coupling one typically encounters in-phase synchro¬ 
nization, at least for sufficiently small Di, precisely where 
the T FPs are correlated. These arguments suggest that the 
approximately constant phase shift in the ensuing time series 
should imply a strong correlation between the TFPs of indi¬ 
vidual units. In other words, the ordering effect of coupling 
can make an impact on the activation processes of units in a 
fashion similar to what is found for long time asymptotic pro¬ 
cesses, such as synchronization between the unit’s time series. 


V. CONCLUSIONS 

We have analyzed the noise-driven first-pulse emission pro¬ 
cess for a single and two interacting type II excitable units 
where both the fast and the slow variable are influenced by 
stochastic perturbations. Our results concern two main issues; 
(j) determining the MPAPs around which the stochastic ac¬ 
tivation paths are clustered, and (ii) examining in detail the 
effects of two different noise sources on the statistical features 
of the activation process, further demonstrating how the statis¬ 
tics is modified due to linear/nonlinear form of interactions. 
For both issues, we have highlighted the impact of stochastic 
bifurcation, which underlies the transition from stochastically 
stable fixed point to stochastically stable limit cycle. 

Since the study is focused on the process of first-pulse emis¬ 
sion, one of the requirements has been to provide a clear- 
cut distinction between the spiking responses and the small- 
amplitude excitations. This has been achieved by introducing 
an appropriate terminating boundary set, given by the spiking 
branch of the cubic nullcline. Note that our problem setup 
is different from the earlier numerical studies on pulse trig¬ 
gering, which have introduced the terminating boundary as 
a fixed threshold llm as well as the recent study on a 
single FFIN unit 111^ . where the ghost separatrix has been 
used as a terminating boundary within the generalized escape 
problem. The advantages of our approach lie in that the ter¬ 
minating boundary is analytically tractable, while the adopted 
formulation further facilitates an immediate generalization of 
the activation problem from the case of a single unit to dif¬ 
ferent scenarios with two interacting units. The differences in 
the event statistics between the three mentioned problem se¬ 
tups may depend on the system parameters b and e, as well as 
the noise intensities. 

Regarding point (i), it has been established that the topo¬ 
logical features of the MPAPs qualitatively depend on which 


type of noise affects the system dynamics. This has been 
demonstrated by examining the MPAPs obtained for three 
characteristic ratios of external vs internal noise, reflecting the 
scenarios where a particular noise source or both sources are 
present. For the fixed characteristic ratio, the MPAP profiles 
change under increasing noise intensity. The changes become 
apparent as one approaches the noise values that give rise to 
stochastic bifurcation. In case of two coupled units, we have 
shown that the topology of MPAPs is substantially affected 
by the linear/nonlinear form of interactions. For the linear 
couplings, the respective MPAPs of the two units are iden¬ 
tical, whereby the solution lies comparably close to that of 
an uncoupled unit under the analogous parameter set. For the 
nonlinear couplings, the MPAPs of two units may become 
visibly asymmetrical, depending on the noise intensity. While 
discussing the topological features of the MPAPs, we have 
also indicated a somewhat surprising numerical finding that 
the trajectories of the effective set of Hamiltonian equations 
selected according to a given predefined rule may show strik¬ 
ing similarity to the MPAP profiles for scenarios involving 
both a single and two coupled FHN units. This observation 
requires a more elaborate study, and suggests that a system¬ 
atic theory possibly adopting certain elements of the standard 
Hamiltonian approach to escape problems may potentially be 
derived for the problem of first pulse emission. 

Concerning point (ii), the statistical properties of the first- 
pulse emission process have been characterized by the depen¬ 
dencies of the mean TFPs and the associated coefficients 
of variation on Di and 02- Note that the previous work on 
FHN model has focused on statistics of first exit times in 
presence of a single noise source which has typ¬ 

ically led to the well-known Kramers result ll^ or its mod¬ 
ifications fnl] . Compared to our approach, the treatment in 
ifl^ |T^ is simplified in terms of definition of the considered 
quantities and with respect to the terminating boundary con¬ 
ditions. An important novel result is that r and R show nearly 
universal dependence on Di and D 2 for a single unit, as well 
as for two interacting units. In particular, the mean TFPs 
are found to display three characteristic regimes, whereby the 
transition from the domain of large r values to the plateau 
region can qualitatively be attributed to the stochastic bifurca¬ 
tion. We have determined the noise intensities that give rise 
to stochastic bifurcation by introducing the model (fTTT i based 
on a Gaussian approximation for the dynamics of a stochastic 
FHN unit. By carrying out the bifurcation analysis of the 
approximate model, we have obtained the Hopf bifurcation 
curve D 2 {Di) which qualitatively outlines the stochastic bi¬ 
furcation of the exact system. In case of two units, the impact 
of stochastic bifurcation has been found to depend on the form 
of coupling. 

We have further examined the distributions of TFPs over 
an ensemble of different stochastic realizations under fixed 
{Di, 02 ). So far, little has been known about the profile of 
corresponding distributions even for a typical escape problem. 
In fact, the only available analytical result, derived from the 
theory of large fluctuations, indicates that the exit time distri¬ 
bution for the escape problem asymptotically acquires expo¬ 
nential form M- However, the problem of first pulse emis- 
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sion is conceptually different from the typical escape prob¬ 
lem, whereas the noise intensities we consider are small, but 
finite. Still, in case of a single unit, we find exponential dis¬ 
tribution of TFPs under prevailing D 2 , which is consistent 
with the Poissonian process. However, the activation process 
dominated by Di yields a different distribution which shows 
a unimodal, Lorentzian-like profile. For the interacting units, 
the profile of TFP distribution is further influenced by the 
form of coupling. 

Apart from the T FP distribution, the form of coupling af¬ 
fects in a nontrivial fashion the correlation of individual acti¬ 
vation events. As a qualitative explanation, in subsection llV Bl 
we have suggested a link between the correlation of individ¬ 
ual activation events and the synchronization features of the 
time series for the given parameter set. In this context, the 
future research should set the ground for potential application 
of the theory of large fluctuations to the research of stochastic 
synchronization in excitable systems. 

This study has shown how the analysis of activation process 


may be extended from a single excitable unit to different in¬ 
stances of two interacting units. In the forthcoming paper, our 
goal will be to examine the noise-driven first pulse emission 
process in an assembly of excitable units, focusing on whether 
such an assembly may be considered a macroscopic excitable 
element, and if so, what may be the analogies and differences 
compared to the excitable behavior of a single unit. 
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